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Abstract. The collapse of cavities under shock is a key problem in various fields 
ranging from erosion of material, ignition of explosive, to sonoluminescence, etc. 
We study such processes using the material-point-method developed recently in the 
field of solid physics. The main points of the research include the relations between 
symmetry of collapsing and the strength of shock, other coexisting interfaces, as well 
as hydrodynamic and thermal-dynamic behaviors ignored by the pure fluid models. In 
the case with strong shock, we study the procedure of jet creation in the cavity; in the 
case with weak shock, we found that the cavity can not be collapsed completely by the 
shock and the cavity may collapse in a nearly isotropic way. The history of collapsing 
significantly influences the distribution of "hot spots" in the shocked material. The 
change in symmetry of collapsing is investigated. Since we use the Mie-Griineisen 
equation of state and the effects of strain rate are not taken into account, the behavior 
is the same if one magnifies the spatial and temporal scales in the same way. 
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1. Introduction 

The first work on cavitation phenomena dated back to 1894 when Reynolds observed 
such a phenomenon in water flowing through a tube with a local constriction p]. A few 
years later, very destructive actions of cavitation to some industrial systems, such as 
steamship propellers, hydroturbine, etc, were found [21 [HI HI [5]. Then, it was addressed 
and discussed under other various backgrounds [fH [7J [S_l [9], for example, hot spots in 
explosives from which reaction ensues[5], the recent observation of sonoluminescence[10j, 
etc. Now, controlled cavitation erosion has been regarded as a powerful tool for modern 
technologies like steam/ wet laser cleaning [llj. effective salmonella destruction [12], and 
treatment for kidney stonespS]. Among the publications on cavitation in literature, 
most of them resorts to experiments, then theoretical calculations, and only a small 
portion rely on numerical simulations. The physical models used can be put into 
two categories, fluidic ones and solid ones. The numerical tools for the former are 
hydrodynamic codes[14], and for the latter are molecular dynamics(MD) [T5 l [T6 | \T7\ and 
finite element method [IS]. 

Most of the analysis here are based on fluidic models. The most elegant 
mathematical description of bubble motion is one-dimensional [5j. It has long been 
observed that a tongue of liquid may be projected into a bubble under gravity and 
a tongue of metal may be projected into a cavity under shock, which confirm that one 
primary mechanism of erosion by cavitation is mechanical [3l [19]. The observation of such 
an asymmetric collapse stimulates the question as to what kind of boundary conditions 
give rise to it. The first kind of attempt simply considers it to be a geometric flow effect 
of the boundary. The second one regards it as the spallation of one surface towards the 
other by an incident shock, where the upstream face converges towards the stationary 
downstream one in the jet. In either of the two cases, the aspherical collapse implies 
that an impulse is applied, typically by an impact jet, to the surrounding fluid or solid 
in contact [5]. 

From the numerical simulation side, fluidic models on cavitation in metals ignore 
some important characteristics of solid, such as plastic strain, hardening, and properties 
relevant to the deformation history. At the same time, the material under investigation 
is highly distorted during the collapsing of cavities. This causes severe problems in 
computational modeling of such a process. It is well known that the Eulerian description 
is not convenient to tracking interfaces. When the Lagrangian formulation is used to 
describe problems with large displacement and large strain, the original element mesh 
becomes distorted so significantly that the mesh has to be re-zoned to restore proper 
shapes of elements. The state fields of mass density, velocities and stresses must be 
mapped from the distorted mesh to the newly generated one. This mapping procedure 
is not a straightforward task, and introduces errors. To treat with large distortions 
problems, several mixed methods have been developed recently. Among these methods, 
the material-point method (MPM), introduced originally in fluid dynamics by Harlow, 
et al[20j and extended to solid mechanics by Burgess, et al [21J, overcomes the above- 
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mentioned troubles. At each time step, calculations consist of two parts: a Lagrangian 
part and a convective one. Firstly, the computational mesh deforms with the body, and 
is used to determine the strain increment, and the stresses in the sequel. Then, the new 
position of the computational mesh is chosen (particularly, it may be the previous one), 
and the velocity field is mapped from the particles to the mesh nodes. Nodal velocities 
are determined using the equivalence of momentum calculated for the particles and for 
the computational grid. The method not only takes advantages of both the Lagrangian 
and Eulerian algorithms but makes it possible to avoid their drawbacks as well. In 
this work we use a improved MPM|22] to investigate the collapse of cavities in shocked 
metal. 

The following part of the paper is organized as follows. Section 2 describes briefly 
the numerical scheme. Section 3 shows the physical model and section 4 presents 
simulation results. Section 5 makes the conclusion. 



2. The Material-Point Method 



The material-point method [21] is a particle method, where the continuum bodies are 
discretized with N p material particles. Each material particle carries the information 
of position x p , velocity v p , mass m p , density p p , stress tensor cr p , strain tensor e p and 
all other internal state variables necessary for the constitutive model, where p is the 
index of particle. At each time step, the mass and velocities of the material particles are 
mapped onto the background computational mesh. The mapped momentum at node 
i is obtained by m^Vj = ^ p m p v p iVj(xp), where Ni is the element shape function and 
the nodal mass rrii reads m 8 = ^ p m p iVj(xp). Suppose that a computational mesh is 
constructed of eight-node cells for three-dimensional problems, then the shape function 
is defined as 

JVi = i(i + £6)(i + Wi)(i + «i). (i) 

where £,T],s are the natural coordinates of the material particle in the cell along the x-, 
y-, and z-directions, respectively, £i,?7i,<ii take corresponding nodal values ±1. The mass 
of each particle is equal and fixed, so the mass conservation equation, dp/di + pV- v = 0, 
is automatically satisfied. 

For pure mechanical problems the differential equation of balance reads, 

pdv/dt = V-<r + pb, (2) 

where p is the mass density, v the velocity, cr the stress tensor and b the body force. 
The momentum equation is solved on a finite element mesh in a lagrangian frame. The 
weak form of it can be found, based on the standard procedure, 

J n p5v ■ dv/dtdfi + J n 5(vV) • (rdVl - J Tt 5v ■ tdT 

- j n p5w • bdfi = . ^' 

Since the continuum bodies is described with the use of a finite set of material particles, 
the mass density can be written as p(x) = J2p=i m p<K x — x p)> where 5 is the Dirac 
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delta function with dimension of the inverse of volume. The substitution of p(x) into 
the weak form of the momentum equation converts the integral to the sums of quantities 
evaluated at the material particles, namely, 

m l dv i /dt=(f,) int + (f j r t , (4) 
where the internal force vector is given by fj mt = —^2p P Tn p cr p -{yNi)/pp, and the 
external force vector reads fj ext = Nib p + ff , where the vector ff is the contacting 
force between two bodies. In present paper, all colliding bodies are composed of the 
same material, and ff is treated with in the same way as the internal force. 

The nodal accelerations are calculated by Eq. (j3J) with an explicit time integrator. 
The critical time step satisfying the stability conditions is the ratio of the smallest cell 
size to the wave speed. Once the motion equations are solved on the cell nodes, the new 
nodal values of acceleration are used to update the velocity of the material particles. 
The strain increment for each material particle is determined using the gradient of nodal 
basis function evaluated at the position of the material particle. The corresponding 
stress increment can be found from the constitutive model. The internal state variables 
can also be completely updated. The computational mesh may be the original one or a 
newly defined one, choose for convenience, for the next time step. More details of the 
algorithm is referred to [22j 123] . 



3. Physical model 

We consider an associative von Mises plasticity model with linear kinematic and isotropic 
hardening |24J . Introducing a linear isotropic elastic relation, the volumetric plastic strain 
is zero, leading to a deviatoric-volumetric decoupling. So, it is convenient to split the 
stress and strain tensors, a and e, as 

cr = s-Pl,P=- l -Tr(<T), (5) 

e = e + ~6I,6 = ±Tr(e), (6) 

where P is the pressure scalar, s the deviatoric stress tensor, and e the deviatoric strain. 
The strain e is generally decomposed as e = e e + e p , where e e and e p are the traceless 
elastic and plastic components, respectively. The material shows a linear elastic response 
until the von Mises yield criterion, 

~N|=<7 y , (7) 

is reached, where ay is the plastic yield stress. The yield cy increases linearly with the 
second invariant of the plastic strain tensor e p , i.e., 

a Y = oy + £ t an ||e p || , (8) 

where cy is the initial yield stress and E tan the tangential module. The deviatoric 
stress s is calculated by 

8 = T iU«, ( 9) 
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Figure 1. (See Figl.jpg))Snapshots of collapse of a single cavity under strong shocks. 
From black to white the gray level in the figure shows the increase of local temperature 
denoted by the plastic work during the deformation procedure. The spatial unit is pan. 
The unit of work is mJ. (a)t=2 ns, (b) t=5 ns. 



where E is the Yang's module and v the Poisson's ratio. Denote the initial material 
density and sound speed by po and Co, respectively. The shock speed U s and the 
particle speed U p after the shock follows a linear relation, U s = Co + \U P , where A 
is a characteristic coefficient of material. The pressure P follows the Mie-Griineisen 
equation of state 



where e is the internal energy, p = p/p — 1, 7 = 7 po/p, C = po c o> D = C(2\ — 1), 
S = C(X — 1)(3A — 1). In this paper the used material is aluminum. The corresponding 
parameters are p = 2700 kg/m 3 , E = 69 Mpa, u = 0.33, ayo = 120 Mpa, £ tan = 384 
MPa, c = 5.35 km/s, A = 1.34, and 70 = 1.96 when the pressure is below 270 GPa [25J. 

Before the numerical experiments, we validated our code by several benchmark 
simulations. The first one involves a cylinder rolling on an inclined rigid plane. The 
second involves the collision of two elastic spheres. The third involves a copper Taylor 
bar impacting to a rigid wall. The fourth is to simulate the process of the collision 
between four identical spheres. All the simulation results agree well with the analytic 
or experimental ones|22j. 

In our numerical experiments the shock wave to the target material is loaded via 
either of two equivalent means, impacting with another body in the vertical direction 
or by a rigid wall at the bottom. The horizontal boundaries are treated as rigid walls. 
The initial configuration of material particles are set to be symmetric about the central 
vertical line. Such a configuration corresponds also to a real system composed of many 
of the simulated ones aligned periodically in the horizontal direction. In this paper we 
focus on the two-dimensional simulations. 

4. Results of numerical experiments 

As the first step, we set a single cavity in the metal material. Due to the boundary 
conditions mentioned above, such a simulation model corresponds also to a very wide 
system with a row of cavities in it. We study various cases where the shock wave changes 
from strong to weak. In the cases with strong shock, we concern the jet creation and 
the distribution of the "hot spots". When the cavity is close to the free surface, we 
check whether or not there are material ejected out of the free surface. In the cases with 
weak shock, we study the effects of cavity size, distance from the cavity center to the 
impacting face, the initial yield stress of the material, etc, to the collapsing procedure. 

Figure 1 shows a snapshot for a case with very high impacting speed. The width of 
the simulated system is 10 pm. Two metal bodies with symmetric configurations collide. 




(10) 
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Figure 2. (See Fig2.jpg)Snapshots of collapse of a single cavity under a strong shock. 
For the pressure contour, from black to white the value increases. The corresponding 
times at Figs.(a)-(f) are 1.2, 1.6, 1.8, 2.0, 2.2, 12 ns, respectively. The spatial unit is 
mm. The unit of pressure is Mpa. 

Figure 3. (See Fig3.jpg)Configurations with local temperature denoted by the plastic 
work during the deformation. The unit of work is mJ. Figs.(a)-(f) here correspond to 
Figs.(a)-(f) in Fig. 2, respectively. 




d (x um) 



Figure 4. Critical impact velocity versus the thickness of the downstream wall of the 
cavity. The symbols are simulation results and the solid line is the fitting result. 

The speed of each is 1500 m/s. That is to say, the relative impacting speed is 3000 m/s. 
Thus, the initial pressure to one body loaded by the shock is about 30 GPa being less 
than the critical value, 270 GPa. The global procedure can be described as follows: (i) 
When the shock wave arrives at the upper-stream wall of the cavity, plastic deformation 
begins to occur. The shock waves at the two sides of the cavity move forwards to the 
free surface. The shock speed at the two sides is larger than the deforming speed of 
the upper-stream wall of the cavity, (ii) The upper-stream wall continue its collapse, a 
configuration with a turned "C" appears. The material projected into the cavity makes 
a jet. "Hot spof'occurs at the tip region of the jet. (See Fig. 1(a).) (iii) The speed of jet 
increases with time. The tip of the jet catches up, then exceeds the shock waves at its 
two sides, (iv) The jetted material impacts the down-stream wall of the cavity, results 
in a pair of vortices rolling in opposite directions. The "hot spots" appear at the centers 
of vortices. (See Fig. 1(b).) 

In the case where the cavity locates near the upper free surface, if the shock wave 
is strong enough, the material projected into the cavity will become a jet and hit the 
downstream wall and break it. This procedure makes another dynamical picture of 
ejection by shock. Such a behavior has been observed in experiments and has some 
potential applications. We show such a dynamical procedure in Fig. 2, where the shock 
wave is loaded by colliding with rigid bottom wall. The simulated system here is 10 jim 
in width. The initial radius of the cavity is 1.5 jim. The distance between the cavity 
center to the upper free surface is 4.5 jim. The impacting speed of the metal body to 
the rigid wall is 1120 m/s. The corresponding times in Figs, (a)-(f) are 1.2 ns, 1.4 ns, 
2.0 ns, 2.4 ns, 4.4 ns, and 12.0 ns, respectively. Fig.2(a) shows a snapshot where the 
shock wave has passed most part of the cavity. The upper-stream part of cavity has 
been substantially deformed and some material have been projected into the cavity. Up 
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Figure 5. (See Fig5.jpg)Snapshots of collapse of a single cavity under a weak shock. 
From black to white the gray level in the figure shows the increase of local pressure. 
The spatial unit is mm. The unit of pressure is Mpa. (a) t=1.0 ns, (b) t=1.6 ns, (c) 
t=2.2 ns, (d) t=3.0 ns, (c) t=5.4 ns, (f) t=16.0 ns. 

Figure 6. (See Fig6.jpg)Configurations with local temperature denoted by the plastic 
work during the deformation. The unit of work is mJ. 

Figure 7. (See Fig7.jpg)Transition of symmetry of collapsing, (a) Asymmetric in 
vertical direction near the impacting face; (b) Nearly symmetric collapse. The gray 
level in the figure corresponds to the plastic work. The spatial unit in the figure is 
mm. The unit of energy is mJ. 

to the time t =1.6 ns, as shown in Fig.2(b), the cavity has been nearly filled with the 
jetted material; the compressive wave is arriving at the free surface from the two sides 
and rarefactive waves will be reflected back. The reflected rarefactive waves decrease the 
pressure in the influenced region and cavitation occurs around the region of the original 
cavity and at the two sides of the jet path. (See Figs. 2(c)-2(d).) Compared with those 
at the two sides, material in the middle is in a much higher pressure and has much 
more pronounced kinetic energy so that it is significantly distorted. The newly created 
cavities coalescence and become a larger one with time. (See Fig.2(e).) If the upper 
wall of the newly created cavity possesses enough kinetic energy, it will be be broken. 
(See Fig. 2(f).) Fig. 3 shows the corresponding configurations with temperature. In 
Fig. (a) the temperature of the compressed region around the cavity is higher that in 
other region. The "hot spot" is at the tip of the metal tongue. In Figs, (b)-(c) the 
"hot spot" occurs at regime hit by the metal tongue. After the appearance of the new 
cavity "hot spots" locate at the inner wall of the cavity, especially the upper and bottom 
walls. (See Figs, (d)-(f).) If there are metal material ejected from the upper free surface 
depends on the initial impacting speed u and the width of the upper wall of the original 
cavity d. The critical impacting speed u increases parabolically with d. Fig. 4 shows 
the simulation results and fitting curve. 

With the decrease of impact speed, the collapse becomes slower. When the impact 
speed decreases to about 200 m/s, the dynamical picture has been significantly different: 
The cavity can not be collapsed completely, and the final configuration around the 
cavity varies from (nearly) symmetric to asymmetric if we change the distances from 
the cavity center to the impacting surface. In Fig. 5 we show an "abnormal" asymmetric 
dynamical picture where the cavity collapsed less in the initial shock direction. Here 
the width of the simulated system is 10 fim. The initial radius of the cavity is 1.5 [im. 
The distance between the cavity center to the impact surface of the metal body is 1.9 
\im. To understand the observed behavior, we "recover" or "magnify" the system in 
such a way: The rigid walls at the two sides and at the bottom of the computational 
region can be regarded as "mirrors". In other words, the system is symmetric about 
the "mirrors". The distance between two successive cavities in the horizontal direction 
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Figure 8. Effects of cavity size [(a) and (b)] and initial yield [(c) and (d)] to the 
collapsing procedure, (a) Area of cavity versus time for different initial radius. What 
shown in the legend are the initial radii, the unit is fi m; (b) Collapsibility versus initial 
radius of cavity; (c) Area of cavity versus time for different initial yield stresses. What 
shown in the legend are the initial yield stresses, the unit is Mpa. (d) Collapsibility 
versus initial yield stress. 



is du = 10/im; while the distance between the cavity in the computed body and the 
fictitious one symmetric about the impacting face is dy = 3.8yum. The cavities reflect 
rarefactive waves. Therefor, the pressure in between decreases. Since the distance dy in 
vertical direction is much less than du in horizontal direction, the rarefactive waves in 
between the two cavities in vertical direction reflected more frequently. Correspondingly, 
the pressure in this region is much lower than those in the surrounding region. This 
results in the less collapsing of the lower part of the cavity. When the shock waves arrive 
at the upper free surface, rarefactive waves are reflected back towards the cavity. This is 
a second reason for the cavity to collapse less in the vertical direction. Compared with 
those from the fictitious cavity below the bottom, the reflected rarefactive waves from 
the upper free surface is much wider. This explains why the collapsing of the lower part 
of the cavity is more pronounced. Fig. 6 shows the corresponding configurations with 
temperature. The hottest region appears in the region below the cavity. In Figs. 5 and 
6, Figs, (a)-(f) correspond to 1.0 ns, 1.62 ns, 2.2 ns, 3.0 ns, 5.4 ns and 16 ns, respectively. 
From Figs. 5(a) and 6(a) it is clear that, around this time, although the rarefactive wave 
decreased the pressure within the influenced region, but the temperature there is higher 
than those in other parts of the system. The reason is that the waves there did more 
plastic work. Since the cavity here corresponds to vacuum, pure fluidic models are not 
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able to capture such a phenomenon. 

The temperature in the "hot spot" increases in the course of collapsing. If we further 
decrease the distance between the cavity and the lower impacting face, the lower part of 
the cavity will be collapsed more pronouncedly. Fig. 7(a) shows the final steady state 
for a case where the lower boundary of the cavity just locates at the impacting face. In 
contrast, if we increase the distance, the collapsed cavity will be more symmetric. Fig. 
7(b) shows a case where the collapsing is nearly isotropic. 

We studied cases with various sizes of cavity. It is found, for the investigated cases, 
that the collapsibility becomes larger when the initial size of the cavity is reduced. The 
variation of area of the cavity with time is shown in Fig. 8(a), where four kinds of radius 
are used. They are 2.5 fim, 2.0 /im, 1.5 /im, and 1 fim, respectively. If we define the 
collapsibility as $ = (Aq — A)/Aq, then it is clear that $ decreases as the cavity becomes 
larger, where Aq and A are the areas of the cavity in the initial and final states. (See 
Fig. 8(b).) In order to study separately the effect of individual material characteristics 
to the collapsing procedure, we changed the initial yield stress of the material under the 
condition that all other parameters are fixed. The corresponding collapsing procedures 
are shown in Fig.8(c), where the numbers in the legend indicates the used initial yield 
stresses. The unit is MPa. Fig. 8(d) shows that the collapsibility decreases parabolically 
if the initial yield becomes larger. 

5. Conclusion 

The collapse of cavities under shock is investigated by the material-point method. What 
mainly concerned in the present paper include behaviors ignored by pure fluidic models 
and those quick distortion procedures which are generally difficult for the traditional 
finite-element method [18J. We studied the relations between symmetry of collapsing 
and the strength of shock, other coexisting interfaces, as well as hydrodynamic and 
thermal-dynamic behaviors. In the case with strong shock, the procedure of jet creation 
in the cavity in studied. The jetted material hits the downstream wall and results in 
two vortices rolling in opposite directions. The "hot spots" occur in the two centers of 
vortices. When the shock is strong enough, the jet material with high kinetic energy 
hits and breaks the downstream wall. The critical impact speed for such a phenomenon 
increases parabolically with the thickness of the downstream wall. In the case with weak 
shock, it is found that the cavity can not be collapsed completely and the cavity may 
collapse in a nearly isotropic way. The transition of symmetry in the course of collapsing 
is relevant to the size of the cavity, the distances to the walls, as well as the impacting 
speed, etc. The distribution of "hot spots" in the shocked material changes for different 
collapsing procedures. Since we use the Mie-Griineisen equation of state and the effects 
of strain rate are not taken into account, the behavior is the same if the spatial and 
temporal scales are magnified in the same way. The results obtained in the paper will 
help to understand better the course of ignition of inhomogeneous explosives, fatigue 
and erosion of metal materials, etc. 
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